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This paper discussed the prediction of oscillatory stability condition of the 
power system using a particle swarm optimization (PSO) technique. 
Indicators namely synchronizing (Ks) and damping (Kd) torque coefficients 
is appointed to justify the angle stability condition in a multi-machine 
system. PSO is proposed and implemented to accelerate the determination of 
angle stability. The proposed algorithm has been confirmed to be more 
accurate with lower computation time compared with evolutionary 
programming (EP) technique. This result also supported with other indicators 
such as eigenvalues determination, damping ratio and least squares method. 
As a result, proposed technique is achievable to determine the oscillatory 
stability problems. 
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1. INTRODUCTION 

With the increase of energy consumption in this age, a study on the stability of the power system 
becomes a necessity, especially the oscillatory stability analysis of power systems. This analysis is used to 
predict electromagnetic swing at low frequencies as a result of undisturbed rotor swing. References [1-10] 
pointed out that the stability of the oscillation in the power system is an important issue. As the power system 
operation changes over time, the stability of the small signal in this power system should be tracked online. 
Selected stability indicators are calculated from the data provided over time to track the system. 
These indicators are updated until a constant value is obtained. In this study, damping torque coefficient Ka 
and synchronizing torque coefficient K; are used as stability indicators. Both K; and Ky values must be 
positive so that the system can be classified as stable [8-10]. 

The least squares (LS) method is a technique commonly used to find the K; and Kg values, 
as reported in [8-10]. However, data update requirements are a major weakness of the LS method. 
In addition, this technique also requires a long computation time. Due to these problems, the LS method 
requires monitoring throughout the duration of the swing. Computational intelligence methods are generally 
used to solve problems in power system stability. Optimisation methods include artificial neural 
network (ANN) [11-13], Genetic Algorithm (GA) [14-16], evolutionary programming (EP) [17-21] and 
Artificial Immune Systems (AIS) [22-24]. ANN is the processing systems that inspired by biological neural 
networks that make up the animal brains. Such systems can learn to perform tasks by considering the 
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examples given. On the other hand, GA is a search technique in which the application is based on 
the combination of natural selection and genetic mechanisms. Major characteristic in GA is the crossover and 
mutation operations which able it to produces high quality solution. The set back of GA is it took too long 
time to converge. Meanwhile, EP and AIS are heuristic population-based search methods that used random 
variation, mutation and selection. Something that distinguishes them is AIS also highlighted a process 
called cloning. This paper discussed a new heuristic approach named Particle Swarm Optimization (PSO) 
technique [25-29]. PSO is influenced by the behaviours of schools of fish and flocks of birds. 
It shows performance beyond EP, AIS and GA methods in searching the optimal solution with faster 
computation time. 

This paper proposes an efficient technique for estimating synchronizing and damping torque 
coefficients in solving oscillatory stability problems. Using this technique, K; and Ky values are estimated 
based on information by three generator responses, namely, the changes in rotor speed (4@(t)), the changes in 
rotor angle (40(t)) and the changes in electromechanical torque (47.(t)). The goal is to minimize the error of 
the estimated coefficients. The IEEE 30-Bus system has been chosen to test the online estimation technique 
for K, and Ka. This paper discussed the oscillatory stability prediction in a multi-machine system using PSO. 
A mathematical model for IEEE 30-Bus system for the angle stability assessment is developed. PSO is 
chosen to optimize the objective function, with J as well as K; and Ka as the control variables. Once the J 
value has been maximized, K, and Ky are analyzed, which verify the stability condition of the system. 
The performance of PSO is then compared with EP and LS. Results obtained from the experiment are then 
verified with the minimum damping ratio (€min) [30-32] and eigenvalues (A). 


2. IMPLEMENTATION OF ANGLE STABILITY ANALYSIS 

The IEEE 30-Bus system model is selected to demonstrate the potential of the proposed technique in 
angle stability assessment for a multi-machine system. Six generators, namely, Generators 1, 2, 5, 8, 11 and 
13 are connected to the buses named Buses 1, 2, 5, 8, 11 and 13, respectively. Reference [9] shows 
the parameters of the system. 


2.1. Philips-heffron model for multi-machine system 

The proposed Phillips-Heffron model for the multi machine system is developed and 
shown in Figure 1. The model is developed on the basis of the single machine of the 
Philips—Heffron model [10]. Kz is the damping torque coefficient, H is the inertia constant, J, and Ky, are the 
time constant and circuit constant of the exciter, respectively. wo is equal to 2afo. Ki~Ko and T3 are constants 
that consist of the function related to the operating real and reactive loading, impedance, electrical torque, 
and the excitation levels in the generator. 


Figure |. Phillips-heffron model for multi-machine system 


2.2. Mathematical modelling of philips-heffron model 
Mathematical modelling can be derived for the proposed Phillips—Heffron model for the multi 
machine system shown in Figure | and is presented in the following mathematical equations: 


Aw; / At = (AT nj — ATo; — Kp;4w;)/2H;,, i = 1,...,m (1) 
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A6;/At = Ww Aw;,i = 1,...,m (2) 
Kao; ) i Kay A0; = Cy -~ 
AE,;/At = , ar ; ATjigp, t= 1,..,.m, j=l,..,m, iF fj 
at ( + DjeiCaj4E,j + Avy /AT 01 J J 
(3) 
—Ks 46; + Djei Ks,ij 45; 
Av,;,/At = , i : T,; — Ave; /AT,;, t= 1,..,m, jJ=1,..,m, t#J 
ef eee + Dizi Ke ij4Eqj / Al et Al J J 
(4) 
Ky Ad; — Dj+i Ki ij 45; ) 
AT,; = a : , |, t=1,..,m, j=l,.,m, iF] (5) 
rm Caer — Dysi ko AE; ; J 
Reference [3] showed the details on (1-5). Can be rewritten into matrix form as follows: 
KSA PBA =H am (6) 
X, = (Aw; Ad; AEG Avgi]",i=1,...,m (7) 
U; = [AT,;], i = 1,...,m (8) 


U; and X; are the input and state signal vectors for i generators, respectively. A; is the system 
parameters’ function with i generators. B; is the disturbance matrix. 


2.3. Synchronizing and damping torque coefficients 
The correlation between the change in estimated electromagnetic torque deviation (47vsi(t)) and 
the changes in rotor angle (46,(t)) and rotor speed (4@,(t)) for the ith generator can be expressed as follows: 


AT.;(t) = K,,A6,(t) + KyAq,(t), i = 1, wey ™M (9) 


where Ky and Kay are K; and Ka for the ith generator, respectively, and m is the number of generators. 
The justification of the stability of a linear system can be performed via the estimation of K, and Ka. 
The positive values of K; and Ka will validate the system as stable. If the system has positive Ks and negative 
Ka, then it is defined to be in the oscillatory instability condition. However, if K; and Ky indicate negative and 
positive values, respectively, then the system is considered to be in the non-oscillatory instability condition. 
In general, the system is unstable if either one of the torque coefficients is negative. 

The stability evaluation of a linear system can be predicted with reference to the K; and Ky values. 
A stable system is guaranteed if the K; and Ky values are positive. If the linear system has positive K; and 
negative Ky, then the system is defined to be in the oscillatory instability condition. The effect of 
the oscillatory instability condition can be detected from the increment of the amplitude oscillations of 
the rotor. Non-oscillatory instability occurs if K; and Ky show negative and positive values, respectively. 
This condition can be verified from the steady increment of rotor angle responses. Detail illustration of 
stable, oscillatory unstable and non-oscillatory unstable conditions can be found in [1]. 


2.4. Eigenvalues and damping ratio 
The scalar parameter of eigenvalues, J can be derived as follows [1]: 


(A—ANg =0 (10) 


Here, the n solutions of 24 (=A), 22, ..., dn ) are the eigenvalues of A. The i eigenvalue can be 
stated as follows: 


Ai = 0; + ja; (11) 
where o; and q; are the real and the imaginary part of the i eigenvalue, respectively. If all value of 4 have 
negative real parts, the linear system is considered stable. The damping ratio (é) for the i eigenvalue is 


defined as the following: 
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f= —G;/J0;? + w;? (12) 


The linear system is certainly in stable condition if all damping ratio have positive value. For simplification 
purposes, only the minimum value of the damping ratio, (Gnin) for the linear system is selected to verify 
the result [30-32]. 


2.5. LS method 
LS technique is often used to obtain the minimum value for the sum of the square of the differences 
between A7.(t) and AT.;(t). The error is defined as [8-10]: 


E(t) = AT, (t) A AT.(t) (13) 


Here, A7.(t) and AT.,(t) are the real and estimated electrical torque, respectively. 

The period of troai must be chosen to estimate the correct value for K, and Ka. The different values of 
tiorat Will result in an inaccurate value for K; and Ky. References [8] stated that the suitable value for tioia: that 
makes K; and Kz constant during the oscillation period is the value of the entire oscillation period. In matrix 
notation, the problem can be described as follows: 


AT,(t) = AT,,(t) + E(t) = Cx + E(t) (14) 
C =[Ad(t) Aw(t)] (15) 
x=[K, Kg]’ (16) 


Here AT.(t) and AT.,(t) are the real and estimated electrical torque, respectively. E(t) is the differences (error) 
between AT.(t) and AT.,(t). Detail calculations can be found in [10]. Although the calculated values are 
accurate, the application of the LS method requires a full oscillation period and takes a long time [8]. 
Therefore, a new indicator is necessary. 


3. OPTIMIZATION TECHNIQUES 

Nowadays, artificial intelligence technology (AI) has been widely used in solving power 
system problems. Evolutionary computation (EC) is one of the AI techniques that promotes logical 
representation approaches. EC is a group of global optimization algorithms that has metaheuristic 
optimization properties. Inspired by biological evolution, among the techniques covered in EC are EP, GA, 
AIS and PSO. EP and PSO are selected as optimization techniques in the present study. 


3.1. PSO 

The PSO algorithm is started with initialization, followed by the update of velocity and position, 
fitness calculation, the best position update and convergence test. Detailed explanations of PSO algorithm 
process are as followed. 


3.1.1. Initialization 

In the initialization process of PSO, the value of synchronizing and damping torque coefficient, K; 
and Kg are generated randomly. In the PSO perspective, K; and Ka are called particles and their values are 
called position (or x). For every position that is created, x;, there is a velocity, v;. In the initialization process, 
the velocity is also randomly created in the range [0, 1]. The random positions are then used to calculate 
the fitness, J. In this initialization process, the i" fitness, J; is set as individual best fitness Jip for i particle. 
For the K; and Ky estimation process, one constraint is identified: the calculated J must be larger than 0.5. 
Initialization process is repeated until total number of initial particles, n is achieved. From these N set of 
particles, the maximum fitness of all particles, Jinax is set as the global best fitness, J,. The position for every 
Jip, Jmax and J, is set respectively as the individual best position p;, position with maximum fitness Pima. and 
global best position pg. 
a. Velocity and position update 

After set of particles are selected in initialization level, all n particles are through a process of 
updating the velocity and position, for every particle. The update process of v; and x; for the i” particle at j” 
iteration is in line with the following equations: 
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VG) = ov,G -N+aiG-D-xG-D}+algoG-D)-xG-D} (17) 
XG) =v,G) +xG - 1) (18) 


Here, c; and cz are acceleration coefficients and « is the inertia weight. 
b. Fitness calculation 
Using the new value of v; and x;, the new fitness, J; are calculated for every new n particles. After 
new fitness is calculated, new value of Jinax and the minimum fitness of all particles, Jmin are selected. 
c. Best position update 
With the new set of position, velocity and fitness for n particles, the update process of individual 
best position, p; and the global best position, pz will be performed if the following conditions are met: 
—  IfJj; is bigger than Jip, select J; as Jip, and select x; as p;. If J; is smaller than or equal with Jj, the value 
of Jip and p; are not changed. 
— If Imax is bigger than Jz, select Jinax aS Jz, and select Pmax AS Pe. Else, if Jinax is smaller than or equal with 
J, the value of J, and pz are not changed. 
d. Convergence test 
Convergence test is attended to regulate the stopping criteria of the optimization process. The search 
process will be terminated if the process has reached the maximum iteration number or the difference 
between the value of Jinax and Jinin is very close to 0. The flow chart which represents the PSO algorithm is 
illustrated in Figure 2(a). 
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(a) PSO (b) EP 
Figure 2. Algorithm flowchart 
3.2. EP 


The Evolutionary Programming (EP) is inspired by the theory of evolution based on natural 
selection. Metaphorically, the breeding of a species will produce offspring with some small variations due to 
mutations. With the competition between offspring and parents in finding the suitability of the environment, 
more suitable members will be chosen next generation. This new generation will reproduce, and this process 
repeats until the suitability between the species and the environment is reached. The overall process of EP 
algorithm is illustrated in Figure 2(b). Detailed explanations of EP algorithm process can be found in [10]. 
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3.3. Objective functions 

In the current study, the objective function formulated is based on the differences of 
the electromagnetic and estimated electromagnetic torques of the i generator, AT.i(t) and AT-si(t), 
respectively, as depicted in (21). This difference or error is estimated for the calculation of K; and Ka for 
every generator in the system. The PSO optimization technique is used to minimize the error with K, and Ka 
being the control variables [10]. 


J, = inv (1 + abs ((47.:(¢) = ATesi(€))/4Ta (©) ) ,i=1..,m (19) 
where m is the number of generators. Hence, the objective function can be defined as follows: 
Maximize (Jj) 


From the optimized J value, a decision can be made to identify the angle stability on the basis of 
the K; and Kz values. 


3.4. Algorithm for angle stability assessment 
The calculation process of Ky; and Ka; for the i” generator is conducted repeatedly to estimate 
successfully the maximum value of J;. The following process is implemented: 
a. Calculate AT.,i(t) using 45i(t), Jo(t), and the estimated torque coefficients using (9). 
b. Evaluate J; using (19). 
c. If Ji is smaller than 1.00, then vary the values of Ksi and Kdi and repeat steps a and b with newly 
generated 45,(t) and 4.i(t) sample data until J; reaches 1.00 or all sample data are used. 
Table | tabulates the parameters used in EP and PSO optimization process. 


Table 1. Parameters of EP and PSO 
Techniques EP PSO 
Parameters B=0.09 cr = C2 =0.9 


4. RESULTS AND DISCUSSION 

The achievement of the PSO technique in estimating K, and Ky are conducted via the IEEE 
30-bus system. Generator data for this system can be found in [7].Three samples of data of rotor angle 46(t), 
electrical torque 47.(t) and rotor speed 4@(t) for all six generators are produced in a MATLAB/Simulink 
environment. Two different values of reactive load at Bus 2 are used to simulate various stability cases. 
The values of the reactive load at Bus 2 are chosen in such a way that two scenarios can be emulated, 
namely, stable and critically unstable conditions as tabulated in Table 2. The three responses, namely angle, 
speed and torque deviations for Case 1 are shown in Figure 3(a), 3(b) and 3(c), respectively. In this case, the 
high damping rate for all responses justifies the system as in stable condition. Overall, all responses are fully 
damped 33 s after the simulation started. 


Table 2. Two different loading conditions 
Case 1 (stable condition) 2 (critically unstable condition) 
Reactive load at Bus 2 35 Mvar 210 Mvar 
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(a) Angle deviation (b) Speed deviation 


(c) Torque deviation 


Figure 3. Responses for stable condition for all generators in case | 


Table 3 tabulates the comparison of K;, Ka, J and number of iterations optimized for six different 
generators from EP, PSO and LS method for Case 1. All the three techniques manage to predict the stability 
conditions correctly indicated by the positive values of torque coefficients K, and Ky for all the 
six generators. In this case, PSO method succeeds to calculate fitness value of 1.000 for all generators. On the 
other hand, EP calculated the highest fitness value of 0.8805 for generator G5. From the iteration perspective, 
PSO and EP are close to each other i.e. between 14~18 iterations. Table 4 shows the comparisons of fitness J 
and iteration result for Case 1. From eigenvalues / perspective, all values are negative, meanwhile the 
minimum damping ratio Cnin give positive value. This confirms that case | is a stable case. 


Table 3. Comparisons of EP, PSO and LS method for case 1 


Gen. Tech. Ks Ka J No. Iter. Gen. _ Tech. Ks Ky J No. Iter. 
EP 4.5152 6.7687 0.8204 17 EP 0.3305 0.6409 0.8077 15 

Gl PSO — 3.8295 7.7154 1.000 14 G4 PSO 0.6419 ~—-0.8375 1.0000 18 
LS 3.4122 7.1221 - - LS 0.5012 0.0819 - - 
EP 0.0402 0.2196 0.8218 15 EP 1.7562 6.0567 0.8885 14 

G2 PSO 0.1063 0.2285 1.0000 14 G5 PSO 1.4818 9.1044 1.0000 14 
LS 0.0903 0.0131 - - LS 1.1976 7.4533 - - 
EP 9.1289 11.3642 0.8113 18 EP 4.3506 3.7928 0.8315 16 

G3 PSO 7.6795 9.5289 1.0000 15 G6 PSO 3.7917 ~— 8.6296 1.0000 16 
LS 6.2812 10.3265 - - LS 4.0129 7.7373 - - 


Table 4. The results of A and Emin for case 1 
Emin A 
-25.32774j82.5561, -25.19494j67.3925, -25.1835+4j66.2373, - 
25.1689+j64.4594, 
-25.17274j65.1107, -25.17704j64.7891, -0.0321, -0.2579, -0.12404j14.6853, 
-0.12574j15.4743, -0.1142+4j16.0849, -0.1424+4)17.1449, -0.1409+4j 16.6304. 
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Case 2 is unstable case which supported by the oscillation increment of angle, speed and torque 
deviation for all six generators, shown in Figure 4(a), 4(b) and 4(c), respectively. Case 2 seems to damp in 
the first place, but after the simulation achieved in 15 seconds, the oscillation of the responses increasing 
dramatically until the simulation end. If analyzing the angle, speed or torque responses in short time i.e. 
10-15 s, this case can be considered as a stable case. Unfortunately, if analyzing these three responses in 50 s 
and above, the oscillations seems increasing gradually and obviously the damping will not stop. Based on 
these results, Case 2.2B is considered as unstable case. 


- Tue (a) 


(a) Angle deviation (b) Speed deviation 


t a) 4 ‘ “ ‘ 


Taw te 


(c) Torque deviation 


Figure 4. Responses for stable condition for all generators in case 2 


The value of Ks, Ka, J and number of iterations for Case 2 are shown in Table 5. EP and PSO 
methods managed to calculate the instability conditions for Case 2 indicated by the negative results of Ky for 
all the six generators. On the other hand, LS failed to deliver correct results as this method calculated positive 
values for both K; and Ka for generators G1 and G5. Based on this result, all the three optimization methods 
are highly recommended to assess the stability condition compare to LS technique. Results showed that PSO 
scored perfect 1.000 in fitness value at the end of the simulation process for all the six generators. EP never 
achieved value of 0.9 for this case. In terms of iteration number, PSO become the fastest method, followed by 
EP. From these results it can be said that PSO is the most capable technique in optimizing the highest quality 
of fitness in the calculation process with admissible computation time. 


Table 5. Comparisons of EP, PSO and LS method for case 2 


Gen. Tech. Ks Ka J No. Iter. Gen. Tech. Ks Ka J No. Iter. 
EP 0.5182 -0.7726 0.7960 17 EP -1.1044 -2.3948 0.8631 20 
Gl PSO 0.6769 -0.5762 1.0000 15 G4 PSO -0.9400 -1.7141 1.0000 15 
LS 0.0730 0.3147 - - LS -0.9340 -1.0216 - - 
EP 2.1324 -1.9547 0.8245 19 EP 2.1978 -0.5515 0.8583 20 
G2 PSO 2.5784 -1.9975 0.9985 15 G5 PSO 2.5632 -0.7174 1.0000 16 
LS 2.6721 -1.8122 - - LS 2.6123 0.0022 - - 
EP -1.7429 -2.0883 0.7886 20 EP -1.1044 -2.3947 0.8163 19 
G3 PSO -1.5087 -2.1320 0.9970 15 G6 PSO -0.9400 -1.7141 1.0000 15 
LS -1.0877 -1.9331 - - LS -0.9218 -2.0166 - - 
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Table 6. The results of A and nin for case 2 
Simin A 
-25.38944j8 1.5538, -25.23234j67.3359, -25.2296+4j66.5607, -25.2493+j64.4474, 
-0.0013 -25.20444j65.5412, -25.17834j65.0325, -0.0014, -0.2568, 0.0184+j 14.6126, 
-0.1184+4j15.5915, -0.1250+j17.0943, -0.03644j16.9197, -0.1278+j 16.4524. 
5. CONCLUSION 


This study has discussed the effectiveness of PSO technique in the oscillatory stability prediction in 


a multi-machine system. In this study, the IEEE 30-Bus test system has been selected. Although both EP and 
PSO are capable to predict correctly the stability condition of all cases, PSO is more convincing compared to 
EP. Optimization via PSO produces higher accuracy for all cases compared with EP. From the iteration point 
of view, PSO and EP are almost the same. 
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